function [beta_cf] = betacf(a, b, x)
% betacf - continued fraction for incomplete beta function computed
% by modified Lentz's method
% From Press et al., p 227
% coded, P. Manis 12/2/98
%
maxit = 100;
eps = 3.0e-7;
fpmin = 1.0e-30;

qab = a+b;
qap = a+1.0;
qam = a-1.0;
c = 1.0;
d = 1.0-qab*x/qap;
if(abs(d) < fpmin) d = fpmin;
end
d = 1.0/d;
h = d;
for m=1:maxit
   m2 = 2*m;
   aa = m*(b-m)*x/((qam+m2)*(a+m2));
   d = 1.0+aa*d;	% one step (the even one) of the recurrence
   if(abs(d)<fpmin) d = fpmin; end;
   c = 1.0+aa/c;
   if(abs(c) < fpmin) c = fpmin; end;
   d = 1.0/d;
   h = h*d*c;
   aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
   d = 1.0+aa*d;	% next step of recurrence (the odd one)
   if(abs(d) < fpmin) d = fpmin; end;
   c = 1.0+aa/c;
   if(abs(c) < fpmin) c = fpmin; end;
   d = 1.0/d;
   del = d*c;
   h = h * del;
   if(abs(del-1.0) < eps) 
      beta_cf = h;
      return;
   end; % are we done?
   
end
if(m > maxit) dips('a or b toob g or maxit too small in betacf');
end
beta_cf = h;
return;
